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ABSTRACT 

We introduce an efficient and accurate alternative to full hydrodynamic simula¬ 
tions, Hydro-PM (HPM), for the study of the low column density Lyman-alpha forest 
{Nm ^ 10^^cm“^). It consists of a Particle-Mesh (PM) solver, modified to compute, 
in addition to the gravitational potential, an effective potential due to the gas pres¬ 
sure. Such an effective potential can be computed from the density field because of a 
tight correlation between density and pressure in the low density limit (5 ^ 10), which 
can be calculated for any photo-reionization history by a method outlined in Hui & 
Gnedin (1997). Such a correlation exists, in part, because of minimal shock-heating in 
the low density limit. We compare carefully the density and velocity fields as well as 
absorption spectra, computed using HPM versus hydrodynamic simulations, and find 
good agreement. We show that HPM is capable of reproducing measurable quantities, 
such as the column density distribution, computed from full hydrodynamic simula¬ 
tions, to a precision comparable to that of observations. We discuss how, by virtue 
of its speed and accuracy, HPM can enable us to use the Lyman-alpha forest as a 
cosmological probe. 

We also discuss in detail the smoothing of the gas (or baryon) fluctuation relative 
to that of the dark matter on small scales due to finite gas pressure: (1) It is shown 
the conventional wisdom that the linear gas fluctuation is smoothed on the Jeans scale 
is incorrect for general reionization (or reheating) history; the correct linear filtering 
scale is in general smaller than the Jeans scale after reheating, but larger prior to 
it. (2) It is demonstrated further that in the mildly nonlinear regime, a PM solver, 
combined with suitable pre-filtering of the initial conditions, can be used to model the 
low density IGM. But such an approximation is shown to be less accurate than HPM, 
unless a non-uniform pre-filtering scheme is implemented. 

Key words: cosmology: theory — intergalactic medium — quasars: absorption lines 
- methods: numerical - hydrodynamics 


1 INTRODUCTION 


tional advances). Several models were proposed to e xplain 
the L yman-alpha. forest (iBahcall fc Salpeter 1965|; Aron: 


The Hpnsitv intprcralartic mpdirim filling- tbp pnormoii.^^ 19^: Black 1984 [Ostriker fc Ikeuchi 1983|; ' Ikeuchi fc Os- 


space fptwppn cralaYTPt; and t.bpTT- ag-grpg-at.innt; nffprs rna- triker 198fj; Rees 1986; [Tkenclh 1986|; Rees 1988|; iBorid, Sza- 


mologis lts a unique and powerful probe of the high redshift 


iinbfprsp (y 2 — .^) st.ill inaccpssiblp for largp galayy siir. 


veys. 


. Thp intprgalactic mpdiiim fbprpaft.pr TCMf manifp^t^ Anuinos, fc Nonuau 1995|; pjernguist et al. 1996|; kliralda- 


lay, & Silk 1988; McGill 199C; Bi, Borner & Chu 1992). How- 
ever, it was only after several groups (ICeneTalMg^ ; Zha^ 


itself of servationallv in the numerous weak absorption lines _ Escude et al. 1996 ; [Wadsley fc Bond 199 6|; ^hang et al. 1997 ) 


along a line of sight to a distant quasar, the Lyman-alpha 
forest. Up to date, an enormous treasury of observational 
data on the Lyman-alpha forest at a wide range of redshifts 
has been collected (see, for example, Hu et al. 1995, Ln et al. 
1996, Cristiani et al. 1996, Kirkman & Tytler 1997, Kim et 
al. 1997 and D’Odorico et al. 1997 for most recent observa- 


performed cosmological hydrodynamic simulations when it 
became apparent that at least an appreciable fraction of the 
Lyman-alpha forest consists of smooth fluctuations in the 
IGM, which arise naturally under gravitational instability, 
rather than discrete absorbers, as has been believed before. 

Abundance of observational data and the existence of 
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2 Gnedin and Hui 


a compelling theoretical framework (i.e. hierarchical clus¬ 
tering) allows one to make detailed comparisons between 
observations and predictions of a given cosmological model. 
Moreover, one might even attempt to use the observational 
data in a maximum-likehood type analysis to infer cosmo¬ 
logical parameters, either in a model-independent way, or at 
least within a class of models. A recent example towards this 
direction is the use of the observed mean Lyman-alpha opti- 
cal depth to put limits on the baryon content o f the universe 
( Miralda -Escude et al. 1996; Rauch, et al.l996 ; Bi fc David- 
sen 19tQ [Weinberg et al. 1997| ). One might contemplate 
going a step further to use other properties of the Lyman- 
alpha forest as a probe of equally interesting cosmological 
parameters such as the normalization and slope of the power 
spectrum (see, for example, Hui, Gnedin & Zhang 1997), the 
massive neutrino density, the epoch of reionization and so 
on. 

However, while hydrodynamic simulations give us in¬ 
sights into the physical properties of the IGM as well as def¬ 
inite predictions for a given cosmological model (provided, 
numerical resolution and physical modelling are adequate), 
computational expense makes them impractical to use in a 
maximum-likehood type of analysis in which a large range 
of models are considered. 

It is therefore important to ask whether a more effi¬ 
cient, and at the same time sufficiently accurate, approxi¬ 
mate method can be developed in place of full hydrodynamic 
simulations. 

Up to date, two different semi-analytical approxi- 
mations have been used : the lognormal appro x imation 
(Bi, Borner & Ghu 1992; Bi & Davidsen 1997; Gnedin 
& Hui jj9^ ) and the (truncated) Zel ’ dovich appro xima- 
tion (Doroshkevich j fc Shandarin 1977 McGill 199C; Hui 


Gnedin, & Zhang 1997). While both approximations are very 
efficient, they might not be sufficiently accurate. For exam¬ 
ple, while the one-point density distribution function is close 
to lognormal for mildly nonlinear fluctuations, the lognormal 
approxim ation itself fails to reproduce the density field ac¬ 
curately (Coles, Melott, & Shandarin 1993). The truncated 
Zel’dovich approximation is somewhat more accurate, and 
can be used for about a decade in the IGM density, from 
about p/3 to 3p, where p is the average density of the uni¬ 
verse. However, a main drawback of the truncated Zel’dovich 
approximation is the necessity of the smoothing of the ini¬ 
tial density held to minimize the amount of orbit-crossing 
by the time of interest. This inevitably introduces artihcial 
smoothing of small scale structure, which could bias one’s 
predictions, depending on the quantities one is interested in. 
While the Zel’dovich approximation can be successfully ap¬ 
plied to stud y the column density distribut ion of the Lyman- 
alpha forest (Hui, Gnedin, & Zhang 1997), it remains to be 
seen whether it can reproduce the detailed results of a hy¬ 
drodynamic simulation with sufficient accuracy. 

In this paper we present a new approximate method, 
which is based on a standard Particle-Mesh (PM) solver, 
modihed to account for the pressure forces acting on a huid 
element. While the PM solver is significantly slower than, 
say, the Zel’dovich approximation, it is still much faster 
than a full hydrodynamic simulation. In order to develop 
a method that will be accurate to within 15% (the reason 
for this number will be clear by §4), we hrst in §2 describe 
two hydrodynamic simulations that we have run to be used 


as templates against which approximate methods are com¬ 
pared. Then, we start with linear theory to develop some 
intuition first. In §3 we discuss the effect of the gas pres¬ 
sure on the evolution of linear perturbations. The conven¬ 
tional wisdom that linear baryon (or gas) perturbations are 
smoothed on the Jeans scale is shown to be incorrect in 
general, and it is demonstrated that the smoothing scale 
depends on the reionization history of the universe. Armed 
with an understanding of the behavior of linear fluctuations, 
in §4 we compare full hydrodynamic simulations with an 
approximate method based on combining a PM solver with 
the appropriate smoothing of initial conditions (with the 
smoothing scale motivated by the linear analysis), as a way 
of taking into account the physical effect of gas pressure 
(this is different from initial smoothing in the case of the 
truncated Zel’dovich approximation as a way of correcting 
for orbit-crossing). We then conclude that this method is 
not sufficiently accurate and proceed to develop our new 
approximation, which we call Hydro-PM (hereafter HPM) 
in §5. 


The idea of HPM is very simple: one modihes a regular 
PM solver to compute, in addition to the usual gravitational 
potential, an effective potential due to the presence of gas 
pressure. This is possible because there exists a tight cor¬ 
relation between temperature and density (or equivalently, 
between pressure and density) in the low density limit, which 
can be c omputed quite accura tely for any given reionization 
history (Hui & Gnedin 1997). A given density field then 
predicts an effective pressure held as well as a dehnite grav¬ 
itational potential held. The fundamental rationale is that 
for the Lyman-alpha forest of sufficiently low column den¬ 
sity (A^hi ^ 10^^ cm~^), shock-heating is not important (or, 
equivalently, the density huctuations are only mildly non¬ 
linear, S ^ 10), which is the one piece of physics that HPM 
does not incorporate. As we will show, this does not com¬ 
promise our accuracy signihcantly while buying us a great 
increase in efficiency over full hydrodynamic simulations. 

It is appropriate that we mention here two wond erful 
pieces of relate d wor k. Petitjean, Mticket & Kates (1995) and 
Mucket et al. (1996) investigated properties of the Lyman- 
alpha forest using PM simulations, suitably modihed to fol¬ 
low the thermodynamics of baryons. Their approach differ 
from our HPM method in at least two aspects: the baryons 
are approximated as following trajectories of dark matter 
prior to shell-crossing (we include the dynamical ehects of 
pressure on baryons), and shock-heating is modelled in their 
method which enables them to study higher column density 
systems. The reader is referred to the above papers for de¬ 
tails. 

Finally, we conclude in §6 with a brief discussion. A 
word on our notation: the symbol p is used to denote mass 
density as usual, as well as the mass density in units of the 
cosmic mean (i.e. p and 1 + 5 used interchangeably). Which 
meaning is intended should be clear from the context. 


2 HYDRODYNAMIC SIMULATIONS 

We have performed two cosmological hydrodynamic sim¬ 
ulations against which we will measure the performance 
of our approximate methods. We used the SLH-P^M code 
as described by Gnedin (1995), Gnedin (1996), Gnedin & 
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Hydrodynamics of the IGM 3 


Table 1. Cosmological Models 


Model h cell size 

LCDM 035 0.055 0 7 065 079 6.6h-^ kpc 

SCDM 1.0 0.05 0.5 0 0.70 15.6/i-i kpc 



1 +Z 

Figure 1. Evolution of the ionizing intensity J 21 for the LCDM 
model {solid line) and the SCDM model {dashed line) as a func¬ 
tion of redshift. 


Bertschinger (1996) and Gnedin & Ostriker (1997). Table 
^ contains cosmological parameters of the two models. We 
have chosen two different models which enable us to test our 
approximate methods under different conditions. For both 
models we have used 64® baryonic mesh and the softening 
parameter was set to 1/3, which gives us a dynamical range 
of 192. Since we are mainly concerned with modelling the 
low density IGM, 5 ^ 10, we do not need to set the softening 
parameter to a very small value. Our choice of the softening 
parameter, however, does enable full resolution of regions 
with (5 = 10 or lower. The LCDM model is identical to the 


mod el u sed in our Equation of State paper (Hui & Gnedin 
1997 |), i^ xcept for the larger value of the softening parame¬ 
ter, while the SCD M model is close to the model studied by 
Zhang et al. (1997). 

The thermal history for each of our simulations is de¬ 
termined by the evolution of the photo-ionizing background. 
Figure shows the evolution of the ionizing intensity J 21 as 
a function of redshift for both hydrodynamic simulations. 
The redshift evolution of J 21 and the spectrum of radiation 
for the LC DM model was adop ted from the simulation de¬ 
scribed in (Hui & Gnedin 1997). For the SCDM model, we 
have adopted the following form of the redshift evolution of 
J 21 : 


>^21(2) = 2 


1 -I- tanh 1 100 


7-z 
8(1 -b z) 


where J 21 i s defin ed in exactly the same way as Jhi in Hui 
& Gnedin (1997), and we adopt the same spectral shape 


as in Hui & Gnedin (1997), equation (7). This form of the 


redshift evolution of J 21 is close to the sudden reioniza¬ 
tion models discussed in detail in Hui & Gnedin (1997). For 
low 2 , 2 <C 7, the ionizing intensity reaches its asymptotic 
value, J 21 = 1.0, and it drops quickly before the redshift of 
reionization, 2rei = 7. The factor of 100 inside the tanh in¬ 


sures that reionization occurs smoothly in a redshift interval 
A 2/2 ^ 0.01. This transition period is introduced to avoid 
numerical instabilities possible when J 21 increases abruptly 
at the redshift of reionization, as in models of sudden reion¬ 
ization. 

Both simulations have been continued until 2 = 3. 
It is worth pointing out that the simulation box in both 
cases was rather small, lh~^ Mpc for the LCDM model, and 
422h~^ kpc for the SCDM model. At the final redshift, fluc¬ 
tuations at the box size are already nonlinear, and simu¬ 
lation boxes are not representative patches of the universe. 
This fact should have minimal effect on the present work: our 
main goal is to develop an understanding of the relationship 
between the dark matter and the baryons on small scales, to 
help us find an approximation that takes into account both 
gravity and gas pressure in an appropriate manner (on large 
scales, things are simple: dark matter and baryons simply 
trace each other). The key is then to resolve structure on 
the relevant small (mass) scales (as will be explained in the 
next section), rather than having a representative sample of 
the universe on large scales. 


3 LINEAR EVOLUTION OF COSMOLOGICAL 
PERTURBATIONS 

The linear evolution of perturbations in the dark matter 
- baryon fluid is governed by two second order differential 
equations: 

^ + 2H^ = 47vGp{fxSx + hSb), 

^ + 2H^ = 47TGp{fxSx + f,5,) - (1) 

where Sx{t,k) s^nd Si,{t,k) are Fourier components of den¬ 
sity fluctuations in the dark matter and baryons (we equate 
baryons with the cosmic gas in this paper) respectively, 
which have respective mass fractions fx and /s, H{t) is the 
Hubble constant, a{t) is the cosmological scale factor, p{t) is 
the average mass density of the universe, cs{t) is the sound 
speed in the cosmic gas (where the sound speed is simply 
defined by c| = dP/dp, assuming an equation of state that 
relates the P and p), k is the comoving wavenumber and t 
is the proper time. 

In the limit where the dark matter is gravitationally 
dominant, /;, = 0 in equation (^), th e growth of J x is de¬ 
scribed by the familiar factor D+{t) ( Peebles 198C ), which 
coincides with a{t) if the matter density of the universe is 
critical. 

The right hand side of the equation for 5b contains two 
terms: the gravitational force and the pressure force. On 
large scales, in the limit k 0, the pressure force can be 
neglected, and the baryon fluctuation obeys the same equa¬ 
tion as the dark matter fluctuation. Assuming that Sb = Sx 
initially, we have 

Sb{t,k 0) — 5x{t,k 0) oi D+ (t). 

On small scales, k —> 00 , the pressure force is dominant, 
and one would expect that the baryon fluctuation is sup¬ 
pressed compared to the dark matter fluctuation. A charac¬ 
teristic scale, where the two forces are equal, is called the 
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4 Gnedin and Hui 


Jeans scale. We denote the wavenumber corresponding to 
the Jeans scale as kj, 

kj = —^A-KGp. ( 2 ) 

cs 

The Jeans scale is in general a function of time, but for the 
special case when the gas temperature T evolves with time 
as T oc 1/a, the Jeans scale is constant in time. In this case, 
and assuming fb = 0 (i.e. the baryons are gravitationally 
subdominant), so that 5x oc D+, equation (|^ can be solved 
analytically: 


5b{t, k) 


5x {t,k) 

1 + k^/ky 


(3) 


Thus, at small scales, k oo, the baryon fluctuation is 
suppressed relative to the dark matter fluctuation by a factor 
of Note that this assumes effectively a very special 

boundary condition: T oc 1/a at all times. 

Let us now consider a more realistic case. At sufficiently 
high redshifts, before the Compton heating of the baryon 
gas by the CMB becomes inefficient, the evolution of the 
baryon temperature is well described by the T oc 1/a law. 
At redshifts of about 100 (depending on the baryon density 
of the universe), Compton scattering becomes inefficient and 
the temperature of the baryon gas drops adiabatically, T oc 
1/a^. By low redshifts, but before the universe reionizes, 
the gas temperature can be treated as effectively zero (or 
in other words, the Jeans mass associated with the CMB 
temperature is too small to be relevant for the modelling of 
the Lyman-alpha forest). 

The gas temperature rises dramatically after the uni¬ 
verse reionizes, and its subsequent evolution is obviously our 
object of interest. As we will show, equations (|^ and (^ no 
longer provide a correct description of the smoothing and 
time evolution of the gas. 

In order to consider a realistic case, we extract the evo¬ 
lution of the sound speed from our SCDM simulation and 
solve equation (|^ numerically with the obtained form of 

Cs{t). 

Figure 0 shows the linear baryon fluctuation, normal¬ 
ized by the linear dark matter fluctuation, at ^ = 3. Three 
different cases are shown: fb — 0.05 as in the simulation (the 
solid line), /;, = 0 in equation (j^, baryons being treated 
as gravitationally subdominant and therefore their gravita¬ 
tional effect is negligible (the dashed line), and fb = 1, all 
matter being baryonic (the dotted line). The two cases with 
fb = 0.05 and fb=0 can barely be distinguished from each 
other in the figure. The case where all the matter is treated 
as baryonic also gives a very similar result. This point will 
be examined further in §6. 

Let us focus on the open circles for the time being. They 
represent the linear baryon fluctuation as given by equation 
(^), where the filtering scale is the Jeans scale as dehned 
in equation (^ for z = 3. One can see that in this realistic 
case (where T does not evolve as 1/a at all times), Hltering 
of the baryon fluctuation occurs at a smaller scale than the 
Jeans scale, contrary to conventional wisdom. In addition, 
oscillations occur at small scales, and the amplitude of these 
oscillations decay at a rate slower that 1/A:^, contrary to 
equation (|^. 

It is possible to understand this analytically. Let us con¬ 
sider the case where the baryons are gravitationally sub¬ 
dominant, fb = 0. Then the dark matter fluctuation simply 




X 






Figure 2. Comparisons of different filtering. The exact linear 
baryon fluctuation for the SCDM model at 2 = 3 as calculated 
from equation are shown for = 0.05 {solid line), = 0 
{dashed line, almost overlapping with the solid line), and /^ = 1 
{dotted line). Points with different symbols represent different 
filtering of the linear dark matter fluctuation (to approximate 
the linear baryon fluctuation): 1/(1 + k'^/k'^) filtering {open cir¬ 
cles), exp(—Altering {filled circles), l/(l+fc^/fc^) Altering 
{filled triangles), and a hybrid Altering which gives the best fit to 
the envelope of the baryon fluctuation (eq. |20|; stars). 


grows like D+{t) (ignoring the decaying mode). Let us con¬ 
sider expanding the ratio of the baryon fluctuation to dark 
matter fluctuation 5b{t,k)/Sx{t,k) in powers of k^. Retain¬ 
ing only the first two dominant terms in the small k limit, 
and recalling that 5b{t, k = 0) = 6x{t, k = 0), we have: 


Sb(t,k) ^ _ A{t) 2 
5x{t,k) D+{t) ’ 


(4) 


where A{t) is an unknown coefficient to be determined. In¬ 
serting equation (^ into (|^ and ignoring terms of order k"^ 
or higher, we obtain the following equation for A{t): 


+ = |ZJ.(^). 


This equation can be easily solved by: 


(5) 


A{t) = l^ dt'cUt')Dyt') (6) 

where the initial conditions A{t = 0) = dA/dt{t = 0) = 0 
are assumed (i.e. no difference between the baryon and dark 
matter fluctuations at early times). Note that A is pos¬ 
itive, which means the baryon fluctuation is always sup¬ 
pressed, compared to the dark matter, in the low k regime. 
We now introduce the filtering scale, with the corresponding 
wavenumber denoted as kp, by the following expression: 


A{t) 


D+jt) 

klit) ’ 


so that equation 


(|) 


can now be rewritten as 
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Sb{t,k) _ 

Sx{t,k) kf 


(7) 


The following expression for the filtering scale kp can be 
obtained: 


kf{t) 


D+{t) Jo 


kj{t') 


dt" 




( 8 ) 


where we have replaced the sound speed by its expression in 
terms of the Jeans scale at the same moment (eq. [^), and 
dot denotes differentiation with respect to the time t. 

An important conclusion follows from equation (|^. Let 
us rewrite it in the following form, using the median value 
theorem: 


kfit) 


kj{t*) 


D+{t) 


dt' 


{b+{t') + 2H{t')b+{t')) 


dt" 


\t" 


where t, is between 0 and t. The expression in square brack¬ 
ets integrates to 1, and we obtain: 


kpit) = kj{t,), 


(9) 


where t* < t. In other words, the filtering scale at a given 
time is equal to the Jeans scale at some earlier time. In 
particular, if the Jeans scale 1/kj is an increasing function of 
time, which is typically the case for sufficiently low redshifts 
after reionization, the filtering scale l/kp is always smaller 
than th 


Jeans scale . The reverse would be true prior to 
reioniz dtion, no wo will ooo in a moment. - 

The above notion of the filtering scale is, strictly speak¬ 
ing, only applicable in the small k limit, because it is de¬ 
rived based on an expansion in k^ (equation |^). To see how 
well this filtering scale provides a description of the linear 
baryon fluctuation in the high k regime, we show in Fig. ^ 
with filled circles the filtering in the form exp{ — k^/kf) (i.e. 
5b = 5x expl—k^/k%]), where kp is computed from equation 
(^) using the evolution of the sound speed (or in other words 
the Jeans scale) as extracted from the SCDM hydrodynamic 
simulation. 

One can see despite the fact that kp is derived in the 
small k limit, the exponential filtering with kp gives an ex¬ 
cellent fit to the baryon fluctuation even for high k, until 
oscillations take over. We also show with filled triangles the 
filtering of the form 1/(1 + k^/kf) (i.e. 5b = 5x/\i + k^/k^]) 
for the same kp, which gives a worse fit for the high k cut-off 
and, as in the case of exp{ — k^/kf) filtering, does not match 
the envelope of oscillations on small scales. 

Encouraged by the excellent performance of the gaus- 
sian filtering on scale of l/kp in reproducing the exact lin¬ 
ear solution, we now consider a few special cases, where kp 
can be calculated analytically. Let us restrict ourselves to 
an flo = 1 universe, where D+{t) = a{t). For simplicity, 
we will assume that the mean molecular weight of the cos¬ 
mic gas does not change, in which case the sound speed is 
directly proportional to the square root of the gas tempera¬ 
ture. First, we consider the case where the gas temperature 
T is zero before reionization (which occurs at a = Orei), and 
remains constant thereafter: 


y _ f 0, a < arei, and 
\ To, otherwise. 

Computing the integral (pi), we obtain for a > a^. 


1 3 


kf kj 10 


In particular, for a Orei, 


( 10 ) 


( 11 ) 



kp = 


Another instructive example is when the gas tempera¬ 
ture decays as 1/a after reionization. 


T = 


■ 0 , 


a ^ a^, 


and 


.TbOrei/a, otherwise. 

In this case the filtering scale for a > a^. 


is given by 


kf 


1 


l-b2 


ff) 


3/2 


-3- 


( 12 ) 


(13) 


In the limit a a^ei we recover the standard result kp = kj, 
but the asymptote is reached only slowly, and even at 2 = 3 
and for Zrei = 7, we obtain kp = 2.2kj. We emphasize the 
departure of the correct filtering scale from the usual Jeans 
scale is a result of the fact that T above is not assumed 
to evolve as 1/a at all times. The time evolution of T con¬ 
sidered above is partly motivated by reionization models in 
which the originally cool cosmic gas was heated up to a high 
temperature by radiation emitted by sources (stars, quasars, 
etc) that turned on at some high redshift. 

Typically, the gas temperature decays as an in terme¬ 
diate power between a° and a~^ after reionization (Hui & 


Gnedin 1997). We, therefore, conclude that in a realistic case 
one should expect that at 2 ~ 3 the filtering scale of the cos¬ 
mic gas is about a factor of 1.5 — 2.5 smaller than the Jeans 
scale, unless the universe reionized at a very high redshift, 

2rei » 10. 

Another interesting example is the evolution of the 
baryon perturbations before reionization. After recombina¬ 
tion at 2 ^ 1200, the cosmic gas temperature is still cou¬ 
pled to the CBR temperature by Compton heating, and 
therefore evolves as T oc 1/a. At a later time aaec = 
0.01(fl6/i^/0.0125)^^®, Compton heating becomes inefficient, 
and the gas temperature decreases adiabatically, T oc 1/a^. 
Since the Jeans scale decreases with time for an adiabatically 
cooling gas, the filtering scale for the cosmic gas is actually 
larger than the Jeans scale. More precisely, a good approx¬ 
imation to the evolution of the cosmic gas temperature is 
given by the following expression: 


T = 


2.73 K/a, 

2.73 Kadec/a^ 


a < adec, and 
otherwise. 


In this case the filtering scale for a > adec is given by 


kl 


1 


31n(o/adec) -3 + 4 


1 / 2 ' 


(14) 


(15) 


For example, for Qbld = 0.0125, kp = 0.45fej at 2 = 10, and 
in term of masses, the characteristic mass scale on which the 
gas distribution is smoothed, Mp oc l/k%, is about 11 times 
larger than the Jeans scale, Mj oc l/fcj. This result has 
important implications for understanding the formation of 
the first bound objects in the universe. 
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6 Gnedin and Hui 


Next, we turn our attention to the oscillations in the 
high k regime, a behavior we can understand analytically 
for the time evolution specified in equation (^|). We can 
solve equation (^) exactly in this case, assuming once again 
the case of an fio = 1 universe with fb = 0 (i.e. baryons 
being gravitationally subdominant): 


5b{t,k) 1 

r k'^ 

n- 

f “ ) 

5x{t,k) l-hfcVfcJ 

71- — n+ 

0>Tei J 



for a > flrei, where 

5 , 3 /T 
4 4V9 3kj’ 

and where 5x grows as a. Note that since rocl/aata> 
flrei, the Jeans scale kj is constant in time. In the limit 
a S> ttrei and for k sufficiently small, equation © reduces 
to equation (^. 

Let us now consider a fixed final a, and take the large 
k limit. Then both n+ and n_ become complex (but 5b is 
still real), and 5b as a function of k oscillates. However, one 
can see that in the high k limit, the amplitude of these os¬ 
cillations is independent of k. We, therefore, conclude that 
in general 5b/5x has no power-law asymptote in the high k 
limit. It is sometimes claimed in the literature that 5b/Sx 
always approaches an asymptote of k~'^ in the high k limit. 
That statement is only correct if T evolves as a fixed power 
law in a at all times (see Bi et al. 1992 for derivation). The 
simple case above provides an example of departure from 
this property. 

Finally, we emphasize that the two hydrodynamic sim¬ 
ulations described in the previous section have sufficiently 
small cell sizes so that the corresponding correct filtering 
scales (l/kp) are resolved by about 5 mesh cells. This en¬ 
sures that we can meaningfully compare different smoothing 
prescriptions, as explained in the following section. 


4 FILTERING INITIAL CONDITIONS FOR A 
PM SIMULATION 

The linear analysis in the previous section shows that the 
two mass components, the dark matter and the cosmic gas, 
evolve differently on small scales: the dark matter is affected 
by gas pressure only via gravitational interaction with the 
gas, while the gas evolution is directly influenced by the ther¬ 
mal pressure on sufficiently small scales. In order to compute 
this complex interaction in every detail, a two-component 
hydrodynamic simulation is needed. But often the preci¬ 
sion achieved by the full hydrodynamic simulation is not 
required. For instance, current observations of the Lyman- 
alpha forest typically give about 10-15% accuracy for the 
column density distribution. We will attempt to develop an 
approximation, that is significantly faster than a hydrody¬ 
namic simulation, but at the same time gives us results with 
similar accuracy. 

As a step toward this goal, we will concentrate in this 
paper on single-component approximations, i.e. approxima¬ 
tions where the evolution of the cosmic gas is computed 


using only one set of resolution elements (in our case parti¬ 
cles) instead of following both the dark matter and the gas 
separately. This approach is certainly more economical than 
a full hydrodynamic simulation, but the question is: can we 
make it accurate enough? 

It is certainly possible to emulate a gas-dynamic solver 
using a simple dark matter solver in the linear regime. Let 
5^^\t, k) and 5^^\t, k) be the linear solutions to equation (|^ 
for a specific cosmological model. Suppose we are interested 
in the baryon fluctuation at some final moment t = tf (which 
is early enough so that the fluctuation remains linear). Let us 
model the evolution of the baryon perturbation with a dark- 
matter-only solver (e.g. PM), which, in the linear regime, is 
equivalent to solving the first of equations (^) and assuming 
/(, = 0. If we choose the following initial condition for the 
dark-matter-only solver at an early time t = ti'. 

5x{t = U,k) = ^^^5f\t = tj,k), (17) 

it is easy to see that we will reproduce the baryon fluctuation 
in the linear regime a,tt = tf. Since, as we have shown in the 
last section, = tf,k) can be modelled hy 5x{t = tf,k) 

multiplied by some suitable filter, the above initial condition 
is equivalent to smoothing the initial 5x with the same filter. 

One might then hope to model the dynamical evolution 
of the gas by a PM simulation, with the initial conditions ap¬ 
propriately smoothed. In other words, one may try to model 
the gas evolution under the assumption that the gas is in¬ 
fluenced by gravity alone, hoping that the initial filtering 
procedure is sufficient to model the effect of pressure. 

There is no guarantee that this simple-minded method 
would work. After all, our idea of a simple filtering scale is 
derived from linear analysis, while for our applications, we 
are interested in regions of space with overdensity below, 
but reaching up to about 10. In fact, we will show in this 
section that this method works to a certain extent, but is not 
good enough, i.e. it fails to achieve an accuracy of 10 — 15% 
in a point-by-point comparison of density and velocity fields 
against full hydrodynamic simulations. Observationally, in¬ 
teresting quantities such as the column density distribution 
are typically measured with an accuracy of about 10 — 15%. 
As we will show in the next section, this level of accuracy 
requires similar accuracies in the density and velocity fields 
themselves. 

Before we embark upon a quantitative comparison of 
the PM -I- filtering method versus hydrodynamic simulation, 
we have to address one technical point. 

A collisionless (alias “N-body”) numerical simulation, 
such as PM, uses particles to follow the evolution of the 
system. For our applications, it is eventually necessary to 
compute the gas density and velocity as a function of spa¬ 
tial positions. How does one convert a distribution of par¬ 
ticles into, say, the density field? There exist several tech¬ 
niques, but in this paper we will adopt the simplest method 
of assigning the density onto a uniform mesh using parti¬ 
cle weights. Specifically, we will use the Triangular-Shape- 
Clouds (TSC) scheme to assign the particle density onto a 
mesh. This method, however, suffers from numerical noise. 
For example, in a sufficiently underdense region a particle 
might be so remote from its neighbors that the TSC assign¬ 
ment would leave empty regions (zero density) between the 
particle and some of its neighbors. This generates unphys- 
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Figure 3. A scatter plot of the dark matter density vs the gas 
density (in units of respective average densities) in the SCDM 
hydrodynamic simulation at z = 3. Because of finite gas pressure, 
the gas distribution does not reach densities as low as those of 
the dark matter. 


ical structure on small scales. The easiest way to suppress 
this numerically generated structnre is to smooth the resul¬ 
tant density distribution with, for example, a gaussian filter. 
However, since we are trying to achieve an order of 10-15% 
accuracy in reproducing the gas distribution, we ought to en¬ 
sure that we reduce the numerical assignment noise to within 
a couple of percent. In other words, what is the degree of 
smoothing we must apply to the TSC-assigned density dis¬ 
tribution to reduce the noise to, say, 2%? 

In order to answer this question, we have performed 
two PM simulations: a low resolution one with 64^ mesh, 
and a high resolution one with 128® mesh with the same 
number of particles (64®) and identical initial conditions. 
The rms overdensity at the final moment is chosen to be 3 
to allow for development of sufficient nonlinearities. Those 
two simulations therefore should produce the same final den¬ 
sity distribution except that the high resolution simulation 
has twice higher spatial resolution. The two density distri¬ 
butions are then smoothed with a gaussian filter with some 
chosen smoothing scale. By comparing the two simulations 
smoothed with varying smoothing scale, we find that the 
smoothing scale should be at least 3 cells to reduce the nu¬ 
merical assignment noise to within 2 %. Th is conclusion has 
also been confirmed by Bertschinger (1997). 

Therefore, from now on, we will present results of col¬ 
lisionless N-body experiments with the final density and 
velocity fields assigned to the uniform mesh by the TSC 
scheme and then smoothed with a gaussian filter of three 
mesh cells. This procedure is admittedly quite wasteful, 
since it implies that we lose a factor of 1.5 to 3 in spatial res¬ 
olution. The advantage is that it is simple to implement. We 
defer developing a more efficient density assignment scheme 
to future work. 


Figure 4. The average (thin lines) and rms {bold lines) frac¬ 
tional errors for the density fields in PM + filtering simulations 
as compared to the gas density field in a full hydrodynamic sim¬ 
ulation for our SCDM model (see Table 1). The different ap¬ 
proximations are: PM -|- exp{—k'^/kf) smoothing (dotted lines), 
PM + 1/(1 + Jkf) smoothing (long-dashed lines), and PM + 
best-fit smoothing (eq. short-dashed lines). Also shown is a 
comparison between the dark matter density from the hydrody¬ 
namic simulation and the density from a PM simulation without 
any filtering (solid lines). The difference in the Green functions 
in the PM and hydrodynamic simulations induces an about 5% 
error. 


Before we move on to testing various forms of PM -|- 
initial filtering, it is interesting to address the question of 
whether we need any initial smoothing at all, i.e. how much 
the dark matter and the gas densities differ in a hydrody¬ 
namic simulation. Figured shows the scatter plot of the dark 
matter versus gas density for the SCDM hydrodynamic sim¬ 
ulation at 2 = 3. One can see that the difference is signifi¬ 
cant, with the dark matter density being a factor of 3 lower 
than the gas density in the lowest density regions. Hence, a 
pure PM simulation, with no modifications to mimic the dy¬ 
namical effects of pressure, would fail to reproduce the gas 
distribution of the low density IGM with sufficient accuracy. 

We now turn to testing the method of combining PM 
with the filtering of the initial conditions, as stated at the 
beginning of this section. A hydrodynamic simulation is run 
for the SCDM model as described in Table All PM sim¬ 
ulations are performed with 64® particles on a 192® mesh 
for the same model. The choice of the mesh size of 192® is a 
natural one given that the hydrodynamic simulation has 64® 
moving mesh and the softening parameter is set to 1/3 (we 
have in fact tested different mesh sizes, from 64® to 256®, 
and found that the 192® mesh gives, as could be expected, 
much better agreement with the hydrodynamic simulation). 
The pre-filtered initial conditions of the PM simulations are 
chosen to be exactly the same as those in the hydrodynamic 
simulation. 

Figure ^ summarizes our findings. Before we proceed 
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further, a few words are in order on our way of present¬ 
ing comparisons between two three-dimensional fields (say, 
density or velocity fields). The easiest way for such a com¬ 
parison is a scatter plot, similar to one presented in Fig. 
^ However, while a scatter plot is sufficiently illustrative, it 
fails to give explicit quantitative information. We, therefore, 
use the following method of comparing two fields hereafter 
in this paper. For definiteness, suppose we are interested in 
the field Q{x) (which could be density, velocity or the spec¬ 
trum; in the case of spectrum, x would be replace by A the 
wavelength). We denote by Qbxact the field taken from one 
of the two hydrodynamic simulations, and by Qapprox the 
field taken from the approximate computation under consid¬ 
eration. Then we identify all spatial points in the relevant 
hydrodynamic simulation which have the value of Qbxact 
within ±0.05dex from some chosen value Qo, and compute 
the following average: 

[Qapprox — QexactJavg = 

(QAPPROx(a;) - QEXACT(a;))|Qj^^^^.j,(,„)^Qp (18) 

and the rms deviation: 

[Qapprox — QexactJpims = 

^((Qapprox(x) - Qexact(®))")|q^^^^^(^,^q^ (19) 

over the ensemble of Qapprox’s at the corresponding spatial 
points in the approximate calculation. The above quantities 
would then be plotted as a function of Qo - In the case of den¬ 
sity, we use Q = Inp where p is measured in units of the cos¬ 
mic average; for the spectrum, we use Q = F/{1 — Fbxact) 
where F is the transmission; for velocity v, we use Q = v/vm 
where Vm is defined in Figure ^ In all cases, the average 
and rms deviations defined above provide quantitative mea¬ 
sures of the fractional error in the corresponding approxi¬ 
mate method, compared against the hydrodynamic simula¬ 
tion. 

Given the two deviations, the average one, and the rms 
one, which one is more important? The average deviation 
can be interpreted as a systematic error: it measures how 
much the “approximate” density field systematically devi¬ 
ates from the “exact” density field. Obviously, it is desirable 
to reduce this error as much as possible. The rms devia¬ 
tion is more like a random error, and while it is also desired 
to be as small as possible, a larger value of the random er¬ 
ror can perhaps be tolerated. In comparing simulations with 
observations, usually a statistical quantity is computed by 
averaging the results of simulations in some fashion. This 
averaging will reduce the random (rms) error, but may not 
reduce the systematic (average) error. Therefore, as we are 
proceeding with our tests, we will try to reduce the average 
error to about 5%, and then try to reduce the rms error 
as much as possible while keeping the average error small. 
We again emphasize that we will concentrate on the density 
range p ^ 10, and will ignore all possible error induced in 
the high density regions. 

Our ultimate object of interest is of course a comparison 
of the gas distributions between two methods, but let us take 
a look at the dark matter distributions first. There are a few 
interesting observations. 

The agreement between the dark matter density from 
the hydrodynamic simulation and the density from a PM 


simulation with identical initial conditions (no pre-filtering 
for this PM simulation; shown with solid lines in Fig. ^ is 
better than 4 percent on average, and about 5% rms, getting 
to 10% at overdensities about 10. At higher overdensities 
the agreement gets worse. What is the reason the two dark 
matter distributions do not agree, in spite of the fact that 
the formal resolutions of two simulations are matched? This 
disagreement is caused by the difference in Green functions 
used to compute the gravity force in two simufations. White 
the hydrodynamic simuiation has the Green function cor¬ 
responding to the Plummer softening, the PM simulations 
have the Green function corresponding to our specific choice 
of density assignment on the PM mesh. This difference in 
Green functions, which is pureiy methodofogicaf, induces er¬ 
ror of up to 10% rms even for <5 ~ 10 (in particular, stronger 
deviation at higher density is due to the fact that the Plum¬ 
mer Green function is slightly softer than our PM Green 
function). 

Finally, we turn our attention to a comparison of gas 
density distributions. We first consider the simplest vari¬ 
ant of the PM + filtering method: we smooth the ini¬ 
tial conditions with the exp{ — k‘^/k%) filter (i.e. S(k) —> 
5{k) exp[—k‘^/kp]), where kp is given W equation (H) with 
kj related to the sound speed through (0), and the evolution 
of the sound speed simply taken from the hydrodynamic sim¬ 
ulation. Recall that this particular choice of filtering gives an 
excellent fit to the exact finear baryon fiuctuation on iarge 
scales (filled circles in in Fig. |^. One might hope that the 
same form of filtering -|- PM gives an adequate approxima¬ 
tion even in the mildly nonlinear regime. 

This case is shown by the dotted lines in Fig. ^ Note 
that while the average error is small for 0.5 p ^ 10, it 
gets significantly worse at p ~ 0.1, and the rms error is 
as high as 20% almost everywhere. What causes the strong 
differences in iow density regions? More specificaiiy, in such 
regions, why does the hydrodynamic simuiation predict gas 
densities substantialfy fewer than the PM -|- smoothing ap¬ 
proximation? One possibie expfanation is that the choice of 
initial filtering is incorrect: the exp{—k‘^/kp) filtering under¬ 
estimates the amount of power at high k. From the linear 
analysis shown in Fig. ^ it can be seen that this filter fails 
to take into account extra power due to oscillations in the 
large k limit. 

We therefore try two other variants of the PM -|- filter¬ 
ing method. One is using the 1/(1 + fc^/fcf-llilter (shown as 
filled triangles in the linear analysis of Fig. H). Its results, as 
compared against the hydrodynamic simulation, are shown 
with the long-dashed lines in Fig. This choice of filtering 
gives a slower cut-off at high k compared to the gaussian fil¬ 
ter. The average agreement at low densities significantly im¬ 
proves with this form of filtering, at the expense of, however, 
increased rms error and average error at higher densities. 

Finally, the short-dashed lines in Fig. ^ show the re¬ 
sults of the PM + filtering method with the following filter 
function: 


Mt,k) = - 


-fcVfcl 


(1 + Fk"^ 


( 20 ) 


This choice of filtering gives a very good fit to the envelope 
of oscillations at high k in the linear fluctuations (the star 
symbols in Fig. ^). One can see that the average error still 
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reaches 10 % within the range p < 10 , and the rms error is 
as high as 30% for intermediate densities. 

We have tried quite a few other forms of filtering the 
initial conditions, each giving effectively different amounts 
of power at high k, but none of them reduces the average 
nor the rms error substantially. 

We believe the fundamental flaw of the above PM + 
filtering procedure is that a single uniform smoothing scale 
is assumed for the whole density field. This is adequate in 
the linear regime where spatial fluctuations of the temper¬ 
ature can be ignored in computing the filtering scale (i.e. 
these fluctuations contribute to terms of higher order than 
those in equation 0). But in the mildly nonlinear regime, 
one can no longer ignore such fluctuations. In fact, places 
with hi gher d ensity tends to have higher temperature (Hui & 
Gnedin 1997 ), and hence higher pressure and more smooth¬ 


ing. One then expects the lower density regions, because of 
their lower thermal pressure, to be less smoothed compared 
to the higher density regions (but confined to p 10). A 
uniform smoothing procedure would tend to overestimate 
the density in the lowest density regions. Note that the PM 
part of our procedure does effectively introduce non-uniform 
smoothing, but it does not do so in a way that mimics the 
action of thermal pressure correctly. 

We are not aware of a computationally efficient way 
of performing the necessary variable smoothing on a large 
mesh. Should such an algorithm be invented, the case for 
the PM -b initial-filtering may be reconsidered, but at the 
moment we must admit that this simple method fails to give 
us the desired accuracy in reproducing results of the full 
cosmological hydrodynamic simulation, and we must search 
for something better. 


5 HYDRO-PM APPROXIMATION FOR THE 
COSMIC GAS DISTRIBUTION 


We have repeatedly emphasized in this paper that dynam¬ 
ically, the main difference between dark matter and gas is 
that the latter is subject to thermal pressure on top of grav¬ 
ity. A hydrodynamic code is designed to compute this ther¬ 
mal pressure and in general, there is no other alternative. 
However, in case of the low density IGM, a very useful fact 
can be exploited to our advantage: there exists a tight corre¬ 
lation (to better than 10 %) between gas density and temper- 
ature (and hence pres sure as well) in the low density regime 
(Hui & Gnedin 1997), where shock-heating is not impor¬ 


tant. The density-temperature relation is well-described by 
a power-law equation of state: 


T^Top-'-\ 


( 21 ) 


where To is a constant of the order of 10 "^K, and 7 is typically 
about 1.4 — 1.6. Both To and 7 evolve with time in a way 
that depends on reionization history, but we have develo ped 
an efficient met hod to predict them with high accuracy ( |Hui| 
fc Gnedin 1997 ). 

The equation of state given above immediately provides 
us with the thermal pressure once the gas density is known. 
The need in the hydrodynamic solver suddenly evaporates, 
and the gas evolution can now be followed with a PM-type 
solver, provided it is modified appropriately to include the 
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effect of thermal pressure. We show below how this can be 
done. 

Let us consider the equation of motion for a cosmic gas 
element: 

^-bRu =-V<()-ivP, (22) 

at p 

where v is the gas peculiar velocity, <!> is the gravitational 
potential, and P is the thermal pressure. If the gas is highly 
ionized (so that the mean molecular weight is roughly con¬ 
stant, which is true for the Lyman-alpha forest), and the 
temperature is a function of density only, so that P = P{p), 
equation (|^) can be reduced to the following equation: 

^ -b Ru = -Vi/> (23) 

where 


Ip = cj> + K, 


(24) 


and Ti, called the specific enthalpy, is 


n{p) = ^ + 


Pip') dp' 


Equation (jS^ is identical to the equation of motion for the 
collisionless dark matter except that the usual gravitational 
potential (f is replaced by an effective potential ip, which 
takes into account both gravity and thermal pressure. Since 
the gravitational potential (p has to be computed from the 
density field in a regular PM simulation anyway, it adds only 
a modest computational overhead to compute the enthalpy 
as well. It is extremely simple to modify the regular PM 
routine to do so, and we will call this method the “Hydro- 
PM”, or HPM. 

In principle, one should then follow the motion of two 
sets of particles: the gas which follows the equation of mo¬ 
tion as in (^) and the dark matter which obeys the same 
equation except that R = 0. In practice, to reduce the com¬ 
putational cost, we treat both sets of particles as if they 
all follow the same equation of motion (equation [^, with 
the full Ip including both gravity and pressure). This might 
seem quite unjustified. But one should bear in mind that on 
large scales, pressure is not dynamically important, and so 
allowing pressure to also act on the dark matter particles 
makes practically no difference. The same cannot be said 
for small scales: the artificially imposed pressure on dark 
matter causes its distribution to be less clustered than it 
should be. It then becomes a question of how sensitive the 
small scale pressure (which is dynamically more important 
than gravity on the same scales) on the baryons is to the de¬ 
tailed distribution of matter. The answer seems to be: not 
very much, but we would come back to this point in the 
last section. For now, the reader can take this single compo¬ 
nent HPM method as a plausible approximation, the merits 
of which can only be weighed through detailed comparisons 
with hydrodynamic simulations. 

There is however an important technical point that we 
should discuss before going onto tests of the HPM method. 
In a PM (or HPM) code, the density is assigned onto a mesh 
using the TSC assignment scheme, as an intermediate step 
in the computation of the potential (p (or ip). As we pointed 
out at the beginning of the previous section, this induces 
numerical noise on small scales (high k). This noise is not 
significant for the gravity calculation, since it is suppressed 
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Figure 5. The average {thin solid lines) and rms (bold solid 
lines) fractional errors for the density fields in HPM simulations 
as compared to the gas density fields in full hydrodynamic sim¬ 
ulations for SCDM and LCDM models, and for different stages 
of evolution, as labeled for each panel. For comparison, two vari¬ 
ants of the PM + filtering method described in §3 are shown: PM 
-I- exp(—smoothing (thin and bold dotted lines for aver¬ 
age and rms deviations compared with hydro) and PM + best-fit 
smoothing {thin and bold dashed lines for average and rms de¬ 
viations compared with hydro! (see Fig. The corresponding 
linear rms overdensity ap (eq. |25|) is also shown for each panel. 


by k~^ power in the computation of the gravitational poten¬ 
tial (0(fc) oc k~^S{k)). The computation of the gas enthalpy, 
however, does not include such suppression, and the numer¬ 
ical noise could be a problem. We therefore smooth the gas 
density according to the prescription (over three mesh cells) 
developed at the beginning of the previous section (in other 
words, we smooth the density field not only at the final 
moment, but also at the intermediate steps of the force cal¬ 
culation). As a result, the pressure force is suppressed on 
scales below about three cell sizes. It is then important that 
we resolve the scale l/kp by at least three cells (assuming 
the linear filtering scale l/kp gives the approximately cor¬ 
rect scale over which the density field is physically smoothed 
due to pressure). Otherwise, the artificially reduced pressure 
at scales below three cells (because of our smoothing pro¬ 
cedure to reduce numerical noise) could lead to unphysical 
clustering on those scales. 

In the limit when the filtering scale is very small, and is 
below the cell size, the pressure effect will be insignificant. 
One then may consider running just a pure PM simulation 
to avoid the additional computational expense of about 25% 
because of the HPM modification. 

Let us proceed to the comparison of the HPM approxi¬ 
mation with full hydrodynamic simulations. We extract the 
equations of state as a function of redshift from our hydro- 
dynamic simulations and use them in the HPM simulation 
(the equations of state thus obtained agree very well with 
those obtained using the method of Hui & Gnedin 1997; we 


use for HPM the exact equations of state from the hydrody¬ 
namic simulation so that we can focus on the error induced 
by the approximate dynamics in HPM). Figure ^ shows the 
average and rms errors for the HPM vs full hydrodynamic 
simulation for the SCDM model at three different epochs 
and for the LCDM model at 2 = 3. We also show for each 
panel the corresponding value of ap, which is the rms linear 
overdensity for the exp{ — k^/kj?) filter: 


2 

ap 


27r^ 


dkk^Ppik, a) exp{ — 2k^/kp{a) 


(25) 


where Ppik, a) is the linear power spectrum of a given model 
at a given value of the scale factor a, l/kp{a) is the filtering 
scale at the same moment given by equation (^) {ap grows 
slower that a because kp increases with time), and the factor 
of 2 in the exponential comes from relating 5 to the power 
spectrum by Ppik) oc 5^{k). The quantity ap therefore mea¬ 
sures the degree of nonlinearity of the gas distribution in the 
model. At 2 = 3, the SCDM model is at a more nonlinear 
state than the LCDM model. 

We also show in Fig. ^two variants of the PM -|- filtering 
methods from Fig. ^ for comparison. 

Note that HPM gives a significantly better fit to the gas 
density distribution than the PM -|- filtering approach. For 
5 ^ 10, the average error generally stays within 5%, and the 
rms error is only weakly dependent on density and is about 
15% for high ap cases, and falling to about 10% for low ap 
cases Q This is an important improvement over the PM -|- 
filtering method. 

The gas density is not the only quantity that we would 
like to model. For the purpose of generating absorption spec¬ 
trum, it is important that we have sufficiently accurate ve¬ 
locities as well. Figure ^ shows the comparison between one¬ 
dimensional gas velocities (velocities projected along some 
fixed direction) in the HPM approximation and in the fnll 
hydrodynamic simulation for our SCDM model (the solid 
line). The quantities on the y-axis in Figure ^ are supposed 
to reflect the average and rms fractional errors in the veloc¬ 
ity. The division by cr„ for small [vexactI is implemented 
to avoid arbitrary blow up of the fractional error when the 
velocity vanishes. The HPM approximation reproduces the 
gas velocity again to within 15% rms error, but the average 
(systematic error) has now increased to more than 10% for 
velocities in excess of two sigma. This is an expected result, 
since high velocities generally correspond to the high density 
regions, where the HPM approximation breaks down (be¬ 
cause shock-heating destroys the tight correlation between 
density and temperature/pressure). We also show for com¬ 
parison results of the PM -|- filtering approximations, which 
cannot quite match the performance of HPM. 

Since we plan to apply the HPM approximation to 
model the Lyman-alpha forest, we must also verify that no 
significant systematic error is introduced in the absorption 
spectra themselves. We generate spectra along randomly ori¬ 
ented lines-of-sight through the hydrodynamic and the HPM 


* The LCDM model shows slightly worse agreement at 5 ^ (5 ^ 
10. This is mostly due to the fact that we saved fewer inter¬ 
mediate data while running this simulation, and as the result, 
the evolution of the equation of state from this simulation is de¬ 
termined less accurately than the respective evolution from the 
SCDM simulation. 
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Figure 6. The average {thin lines) and rms {bold lines) fractional 
velocity errors for: HPM {solid lines), PM + exp(—ini¬ 
tial smoothing {dotted lines) and PM -h best-fit initial smoothing 
(eq. [ 2 ^ dashed lines) as compared against the full hydrodynamic 
simulation for the SCDM model at 2 : = 3. 



0.98 

0.96 

0.94 

0.92 


4863 4864 


\ (A) 



0 0.2 0.4 0.6 


X Mpc) 


Figure 7. A line-of-sight comparison between a full hydrody¬ 
namic simulation {solid line) and the HPM {dotted line) for the 
SCDM model at 2 : = 3. The bottom panel shows the density along 
the line-of-sight, the middle panel shows the peculiar velocity, and 
the upper panel shows the flux as a function of wavelength. This 
line-of-sight goes through an underdense region. 


Figure 8. Another line-of-sight, which goes through an over- 
dense region with <5 < 10. 
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Figure 9. Another line-of-sight, which goes through a highly 
overdense region with 5 > 10. The HPM approximation is ex¬ 
pected to break down in this regime. 


simulations, and show three examples in Figures The 
first line-of-sight passes through an underdense region, the 
second passes through an overdense region with overdensi¬ 
ties S ^ 5 (the HPM method is expected to give accurate 
results in this case), and the third passes through a peak 
with the overdensity 5 ~ 16. The HPM method is expected 
to make a larger error in the third case, and this can be eas¬ 
ily observed in the corresponding bottom and middle panels, 
for density and velocity fields. However, since the transmis- 
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Figure 10. The average {thin lines) and rms (bold lines) frac¬ 
tional decrement errors in an HPM simulation as compared to 
the full hydrodynamic simulations for the SCDM model at z = 3. 
The solid line shows the HPM versus the hydrodynamic simula¬ 
tion for J 21 = 0.3, and the dashed line shows the same comparison 
for J 21 = 0.5. Also shown is the case when the gas temperature is 
decreased by a factor of 100 to reduce thermal smoothing {dotted 
line) when generating the spectra. 


sion F is related to the optical depth t hy F = e“^, and 
the optical depth is in turn approximately proportional to 
density to some power, a relatively large error in density 
prodnces only a relatively small error in F. 

To further quantify this, we show in Figure compar¬ 
isons between the decrements (1 — F) in the full hydrody¬ 
namic simulation and the HPM approximation computed 
from 300 random lines-of-sight. Since the neutral hydrogen 
fraction, and therefore the decrement at a given wavelength, 
depends on the ionizing intensity J 21 , we show two differ¬ 
ent cases: J 21 = 0.3 (the solid line) and J 21 = 0.5 (the 
dashed line). Both values are too high for this model to 
reproduce the observed column density distribution of the 
Lyman-alpha forest. Lower values of the ionizing intensity 
will improve the agreement, since in this case a given value 
of the decrement will correspond to a lower value of the gas 
density. 

One might also wonder if the above comparisons under¬ 
estimate the actual error, because of the small simulation 
box size: the thermal broadening could drastically reduce 
discrepancies, because the broadening width is a fair frac¬ 
tion of the box size in wavelength space. To test this pos¬ 
sibility, we recompute the spectrum for the same lines of 
sight through the HPM and the full hydrodynamic simula¬ 
tions with J 21 = 0.5, but with the gas temperature reduced 
by a factor of 100. The corresponding comparison is plotted 
in Fig. 1^ with the dotted line. One can see that thermal 
broadening cannot explain the small errors in the transmit¬ 
ted flux. 



Figure 11. Column density distributions of the full hydrody¬ 
namic simulation {solid line) and the HPM approximation {dot¬ 
ted line) for the SCDM model at z = 3, computed using the 
Density-Peak Ansatz. A 10% error-bar was added for illustrative 
purpose only. 


Figure ^ clearly shows the range of applicability of the 
HPM approximation. While the average error stays within 
10 %, the rms error is smaller than about 18 % throughout 
the whole range of decrement. A remarkable feature of the 
HPM approximation is that it actually describes regions of 
high decrements rather well. This is because even though 
the HPM method fails to give the right density field with 
sufficient accuracy in high density regions, its errors are ef¬ 
fectively suppressed because the same regions give rise to 
saturated absorption lines. 

As we have emphasized before, in comparing simula¬ 
tions with observations, usually a statistical quantity is com¬ 
puted by averaging the results of simulations in some fash¬ 
ion, which tends to reduce the random (rms) error (in other 
words, the point-by-point comparisons above are a rather 
stringent test). We show one interesting example in Figure 
pd] , namely the column density distribution. We compute the 
column density distributions of the full hydrodynamic simu¬ 
lation and the HPM approximation for our SCDM model at 
z = 3 using the Density-Peak Ansatz (Gnedin & Hui 199f; 
Hui, Gnedin, fc Zhang 1997| ). Both column density distri¬ 
butions are plotted in Figure with the solid and dotted 
lines respectively. We also add a 10% error-bar to the column 
density distribution of the full hydrodynamic simulation for 
illustrative purpose. Note that the two distributions agree 
to within about 13%, and the best-fit slopes differ by less 
than 3%. We thus conclude that the HPM approximation 
can be successfully used to model the Lyman-alpha forest 
when a 10-15% accuracy is sufficient. 
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6 DISCUSSION 

We have demonstrated that the HPM method, based on a 
modihed PM routine to take into account the dynamical 
effect of pressure as well as gravity, is an efficient and accu¬ 
rate alternative to hydrodynamic simulations in predicting 
the density and velocity fields, as well as absorption spectra. 

The key that makes the HPM method possible is the 
fact that in the low density regime, where shock-heating 
is not important, there exists a tight correlation between 
density and temperature/pressure. The almost one-to-one 
relationship between these quantities enables us to rewrite 
the equation of motion of the cosmic gas into a form that 
resembles its collisionless counterpart. The net force on the 
gas is then simply the gradient of an effective potential which 
can be computed from the density field alone. 

The power of the method is enhanced by the fact 
that the density-temperature (or density-pressure) relation, 
which has to be input into the HPM computation, can be 
calculated for any reionization history in a very e fficient 
manne r witho ut running hydrodynamic simulations (Hui & 
Gnedin 1997). 


We have also shown that a somewhat worse accuracy 
(than that of HPM) can be achieved by a simple combination 
of a PM solver and smoothing of initial condition with an 
appropriate filter. 

Both the PM -I- filtering method and HPM use a sin¬ 
gle component model to approximate what is in reality a 
two-component system. The HPM method, in particular, 
treats the dark matter as if it is subject to the same forces 
as the baryons, i.e. gravity as well as thermal pressure. As 
we have explained before, this should not be a problem on 
large scales, because pressure is dynamically subdominant 
on those scales anyway. On small scales, we are indeed in¬ 
troducing an error by allowing pressure to act on the dark 
the dark matter distribution would become less clus — 


matter 

tered titan it slmiild he One fact comes to our rescue how- 
ever: the dominant force on the baryons on small scales is 
thermal pressure, not gravity, and since pressure is deter¬ 
mined by the baryon distribution alone, the actual distri¬ 
bution of baryons on small scales should not be sensitive to 


errors i i the dark matter distribution. The good agreement 
between results of single-component HPM and lull hydro- 
dynamic simulations lends support to this interpretation. 

We can perhaps understand this in a simpler setting. 
In Figu re m, the solid lino shows the baryon fluctuation in 


the limit fb~0 (i.c. when no gravitational effect of the gas 
is included), while the dotted line marks the opposite case, 
fb = 1, (when all the matter is treated as baryonic, or in 
other words, the dark matter is subjected to pressure similar 
to HPM). One might expect quite different behavior between 
the two cases, but in fact they are quite similar. Both on 
large scales and at very small scales (scales of oscillations), 
the fluctuations in both cases almost lie on top of each other. 
It is at the intermediate scales, in fact close to kp, where 
the two depart from each other in a perceptible way. These 
scales however span a rather small range, which is probably 
the reason behind the success of the single-component HPM. 

Finally, a few words on the concept of maximum- 
likehood analysis of the Lyman-alpha forest observations. 
While the HPM method is a factor of 10-100 faster than the 
full hydrodynamic simulations, and only 25% slower than a 


single component PM simulation, it still requires consider¬ 
able computational expense. One can imagine using a more 
efficient, but less accurate, approximation (say, truncated 
Zel’dovich approximation, see Hui et al. 1997) instead of 
HPM. This would introduce larger errors, but will allow us 
to sample a large parameter space of cosmological models. 
When a smaller set of plausible models is crudely identified 
with this technique, one can switch to the HPM and fur¬ 
ther narrow the allowed parameter space to a small region; 
finally, if higher accuracy is desired, several full hydrody¬ 
namic simulations can be run. 

This work was supported in part by the UC Berkeley 
grant 1-443839-07427, and in part by the DOE and by the 
NASA (NAGW-2381) at Fermilab. Simulations were per¬ 
formed on the NCSA Power Challenge Array under the grant 
AST-960015N and on the NCSA 0rigin2000 computer un¬ 
der the grant AST-970006N. 
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